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ABSTRACT 


This  report  describes  the  development  of  a  coordinate  system  generation  pro¬ 
cedure  suitable  for  use  in  the  numerical  prediction  of  the  physical  operation 
of  the  FEATRON.  The  suitability  of  the  coordinate  system  is  demonstrated  by 
application  to  a  simplified  model  of  the  equations  governing  charge  distri¬ 
bution  in  the  vacuum  region  of  a  single  FEATRON  cell.  For  this  purpose,  a 
solution  to  Laplace's  equation  subject  to  mixed  Dirichlet  and  Neumann  boundary 
conditions  is  presented.  The  numerical  method  used  is  the  Linearized  Block 
Implicit  method  of  Briley  and  McDonald  which  is  also  described. 


1.  INTRODUCTION 


The  purpose  of  the  present  study  was  to  develop  a  coordinate  system  generation 
technique  for  subsequent  use  in  the  prediction  of  the  physical  operation  of 
the  FEATRON  shown,  schematically,  in  figure  1.  The  FEATRON  consists  of  (1)  an 
array  of  'cone-shaped*  semiconducting  emitters,  (2)  a  collector  or  anode 
contact,  and  (3)  a  metallic  grid  array  that  controls  the  flow  of  carriers  to 
the  collector  plate.  The  grid  operates  in  a  manner  similar  to  the  grid  in  a 
vacuum  tube,  and  the  gate  in  a  field  effect  transistor.  Transport  within  the 
FEATRON  takes  place  in  the  vacuum,  and  the  device  represents,  through  its 
design,  a  significant  coupling  of  advanced  semiconductor  technology  with  the 
radiation  hard  features  of  the  vacuum  tube.  Indeed,  the  device  itself  must  be 
regarded  as  a  major  advance  in  vacuum  tube  technology  in  which  the  significant 
variation  in  charged  particle  transport  is  associated  with  the  geometry  of  the 
emitter  and  the  placement  of  the  ‘grid-holes'. 

The  present  study  is  confined  to  one  cell  of  the  FEATRON,  shown  schematic¬ 
ally  in  figure  2.  Within  this  cell  transport  is  conceptually  separated  into 
two  distinct  parts:  (1)  carrier  transport  within  the  semiconductor,  and  (2) 
transport  within  the  vacuum  region.  Within  the  semiconductor,  Poisson's 
equation  and  Maxwell's  equations  including  that  of  continuity,  must  be  solved 
along  with  a  set  of  equations  describing  carrier  transport.  Typical  length 
scales  for  the  emitter  are  of  the  order  of  a  micron  or  less  and  nonequilibrium 
effects  offer  the  possibility  of  transient  velocity  overshoot.  Additionally, 
care  is  required  in  treating  transport  associated  with  the  constraints  of  the 
cone  tip,  particularly  with  regard  to  adequately  representing  the  bias 
dependence  of  field  emission.  For  the  vacuum,  transit  is  governed  again  by 
Maxwell's  equations,  coupled  to  momentum,  energy,  and  particle  balance  equa¬ 
tions,  and  relevant  scattering  interactions.  For  the  vacuum,  transit  may  be 
simplified  through  coupling  of  Poisson's  equations  to  kinetic  equations  which 
ignore  scattering. 

The  study  performed  under  Contract  N000 14-83-C-07 19 ,  and  summarized  below  was 
concerned  exclusively  with  the  vacuum  region,  and  represents  the  first  phase 
of  a  multifaceted  study  of  the  charged  particle  dynamics  of  the  FEATRON. 
Specifically,  this  study  was  devoted  to  (1)  developing  a  nonorthogonal  coor¬ 
dinate  system  to  smoothly  represent  the  boundaries  and  interior  of  one  ceil  of 
the  FEATRON,  and  (2)  demonstrating  its  effectiveness  through  solutions  to 
Laplace's  equation  subject  to  mixed  Dirichlet  and  Neumann  boundary  conditions. 

A  clear,  first  extention  of  this  study  would  couple  Poisson's  equation  to 
charge  particle  transport  within  the  vacuum  section  of  the  one  cell  FEATRON. 

This  Report  is  in  seven  sections.  Section  II  contains  a  brief  discussion  of 
the  coordinate  transformation  formalism.  Section  III  is  a  discussion  of  the 
algorithm  used  to  solve  the  problem.  A  description  of  the  geometry  of  the 
boundaries  and  interior  of  the  vacuum  region  is  summarized  in  Section  IV. 
Section  V  displays  the  solution  to  Laplaces  equation,  in  the  absence  of  the 
grid.  Section  VI  introduces  the  modifications  needed  to  include  the  grid. 

The  study  is  summarized  in  Section  VII. 


II.  COORDINATE  TRANSFORMATION 


Laplace's  equation  in  Cartesian  coordinates  is 


A2*  = 


for  the  potential  0  (X,Y,Z).  Due  to  the  assumed  axisymetric  nature  of  the 
FEATRON  cell  simplifications  occur  through  solutions  in  a  two  dimensional 
domain  with  dependent  variables  X[,  X2.  The  choice  of  grid  for  these 
independent  variables  is  critical  for  numerical  calculations.  For  example,  if 
the  independent  variables  X[  and  X2  are  established  on  a  Cartesian  grid, 
the  grid  points  for  the  FEATRON  geometry  will  not  necessarily  lie  on  the 
boundaries.  From  a  numerical  standpoint,  the  potential  distribution  would 
require  interpolation  at  the  boundaries  and  could  incur  spatial  truncation 
errors.  These  numerical  errors  are  avoided  by  using  a  coordinate  system  that 
is  "fitted"  to  the  geometry  of  the  FEATRON.  There  is,  of  course,  a  tradeoff 
in  introducing  curvilinear  coordinates.  The  most  significant  being  that 
equation  (1)  is  more  complicated  in  form. 

The  transformation  of  equation  (1)  to  curvilinear  coordinates  is  direct,  and 
is  found  in  a  variety  of  texts  (see,  e.g.,  [1]).  For  a  set  of  Cartesian 
variables  Xj,  X2,  and  a  set  of  nonorthogonal  curvilinear  variables  Yi, 

Y2  a  set  of  functional  relationships  exist  whereby, 


Not ing 


*t  =  xi  <yt  .  y2> 


*2  =  x2(y,.^ 


dx,  dx9 

=  ^dy>  + 


dx_  <7X, 

d'*  =  ^dy'  *-a^dy* 


a  requirement  on  the  tr ans f ormat ion  of  equations  (1)  and  (2)  is  that  the 
Jacobian,  equation  (6),  does  not  vanish. 


J  = 


-*L  ±L 

<*y,  <*y2 


dx2  dx2 

<*y,  <*y?. 


*  0 


In  terms  of  the  new  variables,  Laplace's  equation  is  written  as 


where 


g1'  = 


tv,  *y, 

d\  dxk 


and  a  summation  convention  is  imposed. 

Equation  (7)  with  a  suitable  set  of  boundary  conditions  is  the  equation  of 
interest  in  this  current  work.  The  numerical  solutions  of  equation  (7)  sub¬ 
ject  to  representative  boundary  conditions  has  been  achieved  using  a  numerical 
algorithm  which  solves  the  steady  state  problem  represented  by  equation  (7), 
by  solving  a  transient  problem  which  reduces  to  the  same  solution  at  steady 
state.  This  solution  procedure  is  discussed  in  detail  next. 


III.  SOLUTION  PROCEDURE 


The  solution  procedure,  summarized  below,  is  based  upon  the  approach  of  Briley 
and  McDonald  [2,3].  Briley  and  McDonald  considered  the  coupled  system  of  % 
nonlinear  time  dependent  multidimensional  equations,  given  by 

=  D(<^)  +  S(^)  (9) 


dH(<^  ) 
~ Jt 


where  represents  a  vector  of  dependent  variables.  For  the  present  study 
4=1,  and  is  the  potential  in  Laplace's  equation.  S  and  H  are  nonlinear 
functions  of  $>£,  and  contain  no  spatial  derivatives.  D($jj,)  is  a  nonlinear 
multi-dimensional  partial  differential  spatial  operator,  as  given  here  by  the 
Laplacian.  For  the  problem  of  interest,  S  is  identically  zero  and  in  the 
steady  state,  the  left  side  of  equation  (9)  is  also  zero,  yielding  a  solution 
to  D ( (f> )  =0 .  To  solve  equation  (9),  H  is  first  backward  differenced  in  time. 
For  a  fully  implicit  scheme 


H"+l  -  H"  =  Dn+'  +  Sn+' 


At 


(10) 


with 


At  =  tn+1  -  tn 


For  a  Crank-Nicolson  formulation 


Hn+*-  Hn  =  D^+D"  +  S^'  +  S" 

At  2  2 


More  generally,  time  differencing  about  tn+^t 


t"  +  /3(tn+  -  tn)  =  /3  tn  +  (l  -  P)\ 


(12) 


n+l  n  n  +  l  n+l  _  ,  n  n 

H  -  H  =  p(D  +  S  )+  (l-/?)(D  +S  ) 


(13) 


At 


The  operator  H,D,  and  S  are  then  formally  linearized  about  tn  as,  e.g., 

in 


d<f> 


n  +  i 


A  +  0(  At^) 


( 14) 


8 


At 


=  P 


°wn  +  s '*>"  *  (fr  +  #) 


dj>  d<f> 


*  (l-/3)(Cl"+s")  U5) 


(w  ~  e7fA')  W"  =  (D“*  s")^  +  0(jf)A+""  Al  <16> 


d<#> 


or 


/dHn 

V  d<#> 


(17) 


Thus  at  each  grid  point,  the  nonlinear  system  of  coupled  PDE's  is  reduced  to 
a  block  of  Ixl  coupled  linear  system  of  temporal  difference  equations,  which 
upon  spatial  differencing  are  solved  once  each  time  step.  Note,  for  the 
solution  to  Laplace's  equation  £=1.  Additionally,  since  the  lineari¬ 
zation  error  is  at  worst  of  the  same  order  as  the  temporal  discretization 
error,  the  linearization  is  not  expected  to  introduce  significant 
inaccurac ies . 

Application  of  equation  (16)  to  second  order  PDE's  in  one  space  dimension 
using  standard  three-point  spatial  difference  approximations  requires  the 
solution  of  one  block  Ixl  tridiagonal  system  per  time  step.  Such  a  system 
can  be  solved  efficiently  using  standard  block-tridiagonal  elimination  pro¬ 
cedures.  However,  application  of  the  LBI  algorithm  given  by  equation  (16) 
through  (18)  to  multidimensional  problems  results  in  the  loss  of  the  narrow, 
block  banded  matrix  structure  obtained  in  one  space  dimension.  The  discre¬ 
tization  of  the  multidimensional  spatial  operator  results  in  a  broad-banded 
structure,  which,  if  solved  by  direct  or  iterative  methods,  can  be  extremely 
inefficient.  Such  observations  led  Briley  and  McDonald  [2,  3]  to  develop 
consistently  split  LBI  algorithms  for  mut id imens  ional  problems.  The  split¬ 
ting  is  accomplished  by  dividing  the  multidimensional  spatial  operator,  L, 


(  18) 


9 


into  one  dimensional  operators  associated  with  each  coordinate  direction. 


L  = 


Lx+  Ly+  Lz 


( 19) 


Equation  (5)  is  then  split  following  the  scalar  ADI  development  of  Douglas 
and  Gunn  (4  J : 


(A+AtLx)A<£*  =  At[D(<£n)  +  S(<£n)] 
(A  +  Atl_y)A<£**  =  AA</>* 

(A  +  AtLz)  A  <£***  =  AA<£** 


In  the  above 


A-  (-£)"- /3A.(-£)" 


d<p 


d<f> 


(20a) 


(20b) 


(20c) 


(21) 


and,  A<J>*,  A<f>**  and  A<)>***  are  intermediate  solutions  of  equations  (20a) 
through  (20c). 

Again,  if  three  point  operators  are  used  to  approximate  the  spatial  opera¬ 
tors,  ,  each  of  equations  (20 a)  through  (20c)  will  be  block 
tridiagonal,  and  can  be  efficiently  solved.  The  block  size  and  band  width 
are  independent  of  the  number  of  grid  points,  hence  the  computational  effort 
required  to  solve  the  sequence  varies  linearly  with  the  total  number  of  grid 
points  regardless  of  the  number  of  space  dimensions  considered.  Elimination 
of  the  intermediate  steps  in  equations  (20a)  through  (20c)  yields 

(  A+AtLK)A~'  (A+AtLy)A-1  (A  +  AtLz)  Ac/)***  =  At[o(<£")  +  S(<£")]  +  0(At2)  (22) 


thus 


A<£n+'  =  A  </>***  +  0(  At3) 


(23) 
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However,  for  two  dimensional  problems,  equation  (20c)  is  omitted,  thus 


A  4>n*'  =  A<£**  +  0(At3) 


(24) 


The  development  given  above  presents  a  brief  outline  of  the  LBI  method  used 
in  the  present  investigation.  A  more  detailed  development,  as  well  as  in- 
depth  discussion  of  LBI  methods,  the  linearization  procedure,  and  related 
topics  may  be  found  in  the  article  by  Briley  and  McDonald  [2]. 

It  should  be  noted  that  in  applying  the  above  procedure  to  the  solution  of 
Laplace's  equation  (£=1)  a  scalar  system  results  and  the  single  element  of 
the  A  matrix  as  defined  by  equation  (21)  is  zero.  Thus,  the  procedure 
outlined  above  cannot  be  directly  applied  to  solve  Laplace's  equation. 
However,  this  problem  can  be  circumvented  by  adding  a  pseudo  time  term  to 
Laplace's  equation  and  marching  in  this  pseudo  time  space  to  a  steady  state. 
The  steady  state  solution  thus  generated  satisfies  Laplace's  equation  and  is 
the  desired  solution.  Additionally,  the  pseudo  time  step  can  be  varied  in  a 
manner  which  accelerates  convergence  to  the  steady  state  resulting  in  an 
efficient  solution  technique. 


For  the  more  general  problem  of  transport  within  the  FEATRON,  the  application 
of  the  LBI  procedure  outlined  above  to  the  relevant,  time  dependent  transport 
equations  coupled  to  Poisson's  equation  is  also  not  a  straight  forward  task. 
The  direct  application  of  LBI  procedures  to  this  system  yields  an  A  matrix  of 
the  form 


A 


l 

XXIX 

I 

X  XIX 


6“  olo 


(25) 


which  is  singular,  due  to  the  absence  of  a  time  derivative  in  Poisson's 
equation.  While  this  does  not  create  problems  if  only  one  dimensional  solu¬ 
tions  are  to  be  computed,  for  multidimensional  solutions  this  singularity 
makes  it  impossible  to  apply  the  ADI  splitting  given  by  equation  (20)  to  this 
system,  and  an  alternative  must  be  sought.  One  alternative  is  to  again 
introduce  artificial  time  derivatives  in  each  equation.  For  example,  in 
Poisson's  equation,  the  artifical  time  derivative  yields 

d<f>  , 

-3T  ♦  V  *  =  s  (26) 


The  physical  time  derivatives  in  the  other  governing  equations  are  then 
treated  as  source  terms,  and  the  artificial  time  derivatives  track  the 
development  of  the  solution  variables  in  an  iteration  space.  Such 
"iterative"  LBI  schemes  can  be  highly  efficient  for  generating  steady 
solutions  to  coupled  systems  of 


1  1 


nonlinear  PDE's  since  the  artificial  time  step  can  be  chosen  to  maximize  the 
convergence  rate.  It  is  noted,  however,  that  for  transient  applications  such 
an  approach  would  require  that  a  computational  effort  of  the  order  of  that 
required  to  generate  a  steady  state  solution  be  expended  at  each  physical 
time  step.  This  would  obviously  result  in  highly  excessive  computation  time, 
and  would  be  impractical.  Another  approach  [5],  which  has  been  used  success¬ 
fully,  is  to  reformulate  the  system  of  governing  equations  in  such  a  manner 
that  Poisson's  equation  can  be  decoupled  as  is  suggested  by  the  partitioning 
of  the  A  matrix  shown  in  equation  (25).  This  allows  the  time  dependent 
governing  equations  to  be  solved  efficiently,  without  iteration,  using  an  LBI 
scheme,  followed  by  solution  of  Poisson's  using  the  approach  discussed 
previously  for  Laplace's  equation.  Such  an  approach  can  be  used  to  examine 
transport  with  the  FEATRON,  where  additional  work  should  concentrate  on 
implementing  the  moments  of  the  Boltzmann  transport  equation  for  studying 
carrier  transport  within  the  structure. 
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II 


FIGURE  4,  GRID  DISTRIBUTION  FOR  A  FEATRON  CELL 


BOUNDARY 


SCHEMATIC  OF  THE  FEATRON  DEVICE  SHOWING  THE  RELEVANT  DIMENSIONS  OF  THE  PROBLEM 


FIGURE  2.  SCHEMATIC  OF  ONE  CELL  OF  THE  FEATRON  TUBE 
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VII.  SUMMARY  AND  PROPOSED  FUTURE  DIRECTIONS 


The  objective  of  the  present  study  was  to  design  and  implement  a  coordinate 
system  and  grid  structure  to  represent  the  constrained  design  of  one  cell  of 
the  seraiconductor/vacuum  FEATRON  device.  The  task  was  undertaken  to  enable  a 
numerical  study  of  the  space  and  time  dependent  operational  properties  of  the 
FEATRON.  The  objectives  of  the  study  were  successful. 

To  undertake  future  studies  of  the  transient  characteristics  of  the  FEATRON, 
it  is  necessary  to  establish  a  set  of  relevant  governing  transport  equations. 
In  this  case,  there  are  two  regions  of  study,  the  semiconductor  tip  and  the 
vacuum  region.  It  is  suggested  that  the  next  step  in  the  FEATRON  simulation 
involve  solutions  of  the  transport  equations  within  the  vacuum  region.  Here 
two  approaches  may  be  taken.  The  simplest  is  to  assume  a  kinetic  equation 
within  the  vacuum  that  ignores  scattering,  and  to  solve  this  equation 
s imul t aneously  with  Maxwell's  electromagnetic  equation.  This  approach  has 
been  discussed  in  the  past  and  is  part  of  most  undergraduate  texts.  A  more 
careful  approach  involves  obtaining  solutions  to  the  Boltzmann  transport 
equation  simultaneously  with  Maxwell's  equations.  Workers  at  SRA  have  had  a 
broad  range  of  experience  in  obtaining  approximate  solutions  to  the  Boltzmann 
transport  equation  through  its  moments.  These  equations  include  the  continu¬ 
ity,  momentum,  and  energy  balance  equations  and  incorporate  the  influence  of 
collisons  [6J.  These  moments  equations  are  for  multispecie  transport,  which 
in  the  present  case,  would  incorporate  the  transport  of  emitted  electrons  and 
any  other  charged  part icles -  specif ic  to  the  vacuum  region.  It  ia  suggested 
that  these  equations  be  used  to  study  transport  within  the  FEATRON. 

The  goals  of  such  a  study  would  include  examining  the  time  dependent  large 
signal  behavior  of  the  FEATRON  with  such  output  as  power  and  efficiency.  The 
dependency  of  the  power  and  efficiency  of  the  FEATRON  on  design  parameters, 
for  an  assumed  set  of  boundary  conditions  representative  of  carrier  emission 
from  the  FEATRON  tip,  would  also  be  studied.  The  coordinate  system  to  be  used 
for  this  study  has  already  been  developed  and  would  minimize  numerical  errors 
associated  with  the  highly  constrained  FEATRON  tip. 


SECTION  VI.  ALGORITHM  MODIFICATION  FOR  INCLUSION  OF  THE  METALIC  GRID 


The  only  detailed  calculation  performed  involved  the  solution  to  Laplace's 
equation.  However,  in  anticipation  of  future  studies,  modification  of  the 
grid  generation  strategy  was  introduced  to  incorporate  the  grid.  Figure  6 
shows  an  example  of  a  coordinate  system  in  which  a  FEATRON  grid  has  been 
included.  This  coordinate  system  is  designed  to  concentrate  points  in  regions 
of  anticipated  high  potential  gradient.  The  code  used  to  generate  this  mesh 
has  been  written  to  facilitate  examination  of  alternative  designs  by  allowing 
for  changes  of  the  grid  distribution  through  modification  of  the  code  input 
parameters . 


SECTION  V.  SOLUTIONS  TO  LAPLACE’S  EQUATION 


Laplace's  equation  governing  the  potential  0(x,y)  in  the  field  a  single  cell 
of  the  FEATRON  device  has  been  solved  subject  to  the  following  boundary 
conditions:  (see  also  figure  2) 

(i)  Fixed  potential  on  the  collector  (0“0c) 

(ii)  Fixed  potential  on  the  emitter  (0"0e) 

(iii)  On  all  the  surfaces,  the  first  derivative  of  potential  in  the 
direction  normal  to  the  surface  is  set  to  aero. 


It  should  be  noted  that  in  this  initial  phase  of  the  work,  the  metallic 
control  grid  shown  in  figure  2  has  been  omitted  from  the  calculations.  Also, 
even  though  the  emitter  tip  is  a  semiconductor  and  as  such  would  not  sustain  a 
fixed  potential,  in  the  absence  of  a  suitable  model  for  the  tip  potential 
distribution,  a  fixed  potential  was  assumed.  No  additional  numerical 
difficulties  are  incurred  by  including  appropriate  boundary  conditions, 
obtained  from  a  model  for  the  tip  region  physics. 

Figure  4  shows  the  computational  coordinate  system  used  in  the  current  calcu¬ 
lation.  Figure  5  shows  the  potential  distribution  obtained  by  performing  the 
calculation  described  in  this  report.  In  view  of  the  axisymmetric  nature  of 
the  problem  only  half  of  a  FEATRON  cell  has  been  shown.  In  the  inanediate 
vicinity  of  the  emitter  a  relatively  large  potential  gradient  exists.  The 
field  is  two  dimensional  in  nature  decaying  to  a  linearly  varying  almost  on 
dimensional  distribution  in  the  vicinity  of  the  collector.  This  result  is 
qualitatively  as  anticipated  and  serves  to  illustrate  the  validity  of  the 
approach  developed  under  this  effort. 


The  technique  used  in  this  study  is  to  constrain  the  values  of  the  physical 
coordinate,  y^,  at  specific  values  of  the  computational  coordinate, 


■ck* 


At  interior  points,  the  h^'s  are  referred  to  as  interior 


cluster  points.  In  this  formulation,  two  pivot  points  are  associated  with 
each  interior  cluster  point  and  one  pivot  point  is  associated  with  each  end 
cluster  point  (hpmin  and  npmax^ •  Thus,  if  there  are  k  interior  cluster 
points,  the  total  number  of  pivot  points  is  N=2k+2.  To  determine  the 
functional  relationship  y(q)  described  in  equation  (40),  requires  (if  the 
values  of  npj  and  oj  are  prescribed)  that  the  values  of  the  N+10's  in  that 

equation  be  calculated.  Both  equation  (40)  and  equation  (39)  are  linear 
equations  with  respect  to  the  0's.  By  constraining  the  values  of  the  physical 
coordinates  at  the  interior  and  end  cluster  points,  k+1  linear  independent 
equations  in  the  0's  are  obtained.  By  further  constraining  the  values  of  hte 
slope,  f(n),  at  the  interior  and  end  cluster  points,  a  further  k+2  linear 
indpendent  equations  for  the  N+l=2k+30's.  This  system  is  solved  for  the  0's 
by  standard  Gaussian  elimination  techniques  and,  hence,  all  the  constants  of 
equation  (3)  are  uniquely  determined. 


There  are  several  advantages  to  the  use  of  the  series  of  complimentary  error 
funct ions : 


(1)  f(n)  is  positive  finite  and  non-zero,  i.e.,  y  will  always  increase 
with  increasing  n. 

(2)  f(n)  is  continuous  and  successively  differential  and  integrable;  and 

(3)  If  the  pivot  points  are  spaced  at  a  greater  distance,  (in 
computational  space)  than  a/2  from  each  other,  the  complimentary 
error  functions  will  have  minimal  interaction.  Thus,  the  width  and 
location  of  a  complimentary  error  function  can  be  changed  without 
effecting  other  complimentary  error  functions. 


(c)  Internal  Grid  Distribution 

After  boundary  grid  points  are  distributed,  the  coordinates  within  the 
interior  of  the  computational  domain  are  constructed.  This  is  achieved  by 
connecting  corresponding  points  on  the  top  and  bottom  surfaces  by  straight 
lines.  In  general,  these  line's  are  not  perpendicular  to  the  top  and  bottom 
surfaces,  but  rather  inclined  to  minimize  the  nonorthogonality  at  the  emitter 
surface,  as  seen  in  figure  4. 

The  second  family  of  coordinate  curves  is  then  determined  by  distributing 
along  each  of  the  lines  of  the  first  family  according  to  a  normalized 
distribution  designed  to  concentrate  points  near  the  emitter  boundary  where 
large  potential  gradients  are  expected. 


16  - 


which  can  be  integrated  to  yield 


y  =  /  Hfjldyf  + 


rmin 


(38) 


s(n)  is  chosen  to  be  a  series  of  N  complimentary  error  functions  of  the  form 


HV)  =  ~  =  0O  +  T  I  {erfc[-^-^- v>]  -  ['  +  signtaj)]}  /3§ 

'  j:|  i  1 


(39) 


Ihe  jth  complimentary  error  function  is  centered  in  computational  space  at 
location  rjpj  (which  is  referred  to  as  a  pivot  point),  and  aj  is  the 

width  in  computational  space  in  which  90Z  of  the  grid  size  variation  takes 
place.  At  computational  space  location  np ■ ,  the  value  of  f(n)  assumes  a 
local  maximum.  Y  is  a  scaling  constant  for  a; (Y*l. 63x2«2.326) , 
erfc  (1. 163)=0. 10)  and  in  the  limit  as  aj+0,  8j  is  the  jtl*  step  height 
for  pivot  j,  i.e.,  the  difference  in  grid  spacing  on  either  side  of  qp  ^ . 

Substitution  of  equation  (39)  into  equation  (38),  and  integrating  yields 


y  -  ymin  =  +  I 

-  [l  +  sign 


(40) 


where 


d)W  =  (41) 
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Section  m  (x3  <  x,  <  x4  ;  y3  <  y.  <  y4) 

x(.  =  x3  +  (S.  -  S3)  sin  0 

yi  =  y  3  *  (S.  -  S3)cos  B 

y4  =  H  -  R2(l  -  Sin  B) 

*4  =  x3  +  (y4-  y3) tan  6 

s4  =  s3  +  (y4-  y3)/cos  6 


(32) 

(33) 

[99.79  xIO-6  cm] 
[o.4996  I0~6cm] 


Section  IV  (x4  <  x.  <  x5;  y4  <  y.(  <  y5  ) 

x,  =  x4  +  R2  [cos0  -  co$(0  +<£|)] 

y,  =  y4  +  R2  [sin(0+^>,)  -  sin®] 

0,  =  (S,  -  s4)/r2 

x5  =  x4  +  2R2cos# 

*5  =  U 

S5  =  S4  +  R2(ir  -  20  ) 


(34) 

(35) 

(36) 

[o.5 x  IO-3cm] 

[99.79  x  I0~6cm] 


(b)  Distribution  of  Grid  Points  on  the  Boundaries 

After  the  surfaces  are  analytically  defined,  grid  points  are  positioned  on 
the  boundaries.  The  number  of  grid  points  is  chosen  to  be  sufficient  to 
yield  required  resolution  of  the  potential  and  of  the  boundary  geometry. 
This  requires,  for  example,  clusters  of  grid  points  on  Sections  2  and  6. 


Trans  format  ion  Function 

The  distribution  of  grid  points  on  boundaries  is  accomplished  by  use  of  a 
transformation  technique  described  in  l 7 ) .  If  y  and  n  designate  the 
independent  variables  in  physical  and  computational  (grid  point)  space, 
respectively,  a  transformation  function  is  introduced: 


dy 

dr, 


(37) 
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IV.  DESCRIPTION  OF  THE  BOUNDARIES 


To  achieve  numerically  accurate  and  economical  solutions  it  is  necessary  to 
(1)  describe  the  boundaries  accurately,  (2)  distribute  grid  points  smoothly 
along  the  boundaries,  and  (3)  construct  the  grid  inside  the  computational 
domains.  The  boundary  description  follows. 

(a)  Boundary  Description 

With  reference  to  figure  3,  there  are  four  boundaries  to  consider.  The  most 
complex  is  Boundary  1,  and  is  described  in  detail  below. 

Boundary  1 : 

Section  I  (x,  <  x.  <  x2) 


y,  «  o 

(27) 

*i  =  *.  +  si 

where  Sj  is  the  accumulated  arc  length  along  Boundary  I. 
x,  =  0 

(28) 

Section  II 

(x2  <  Xj  <  x3  ;  y2  <  y.  <  y3  ) 

Vl 

= 

y2  +  R,(l  -  cos  <^j) 

(29) 

*i 

= 

x2  +  R(sin<£(. 

(30) 

4>-, 

r 

<S,  -  S,)/R, 

(3  1) 

where 

is  the  angle  subtended  at  the  center 

of  the  circular  arc 

fillet  by  radial 

lines  through  point  (x2,  y2 )  ond  the  current 

location  (x.  ,  y(  ). 

*2 

- 

x3  -  R,  cos  0 

[42.2058xl0'3cm 

] 

y2 

= 

0 

*3 

= 

xc  -  H  -  ^R2(l  -  sin0)  -  y3jton0 

[43.0l45x  IO-,5cmJ 

R. 

is 

the  radius  of  fillet  (see  figure  3) 

[o.IxlO-6  cm] 

r2 

is 

the  rodius  of  the  tip 

[o.5xlO-6  cm] 

9 

is 

half  the  cone  angle 

[36*] 

H 

is 

the  cone  height 

[lO-4  cm] 

y  3 

= 

y2  +  R,  (1  -  sin0) 

[4  122 lx  I0"6  cm] 

*c 

is 

the  x  coordinate  of  the  o.~e  center  1 

ine 

[0.5  x  I0-3  cm] 
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